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Abstract 

We describe and compare two types of microwave sky simulations which 
are good for small angular scales. The first type uses expansions in spherical 
harmonics, and the second one is based on plane waves and the Fast Fourier 
Transform. The angular power spectrum is extracted from maps correspond- 
ing to both types of simulations, and the resulting spectra are appropriately 
compared. In this way, the features and usefulness of Fourier simulations are 
pointed out. For I > 100, all the simulations lead to similar accuracies; how- 
ever, the CPU cost of Fourier simulations is ~ 10 times smaller than that for 
spherical harmonic simulations. For I < 100, the simulations based on spherical 
harmonics seem to be preferable. 

Subject headings: cosmic microwave background — cosmology: theory — large-scale 
structure of the universe 



1 Introduction 

There are important questions which could be addressed through simulated maps of 
the microwave sky. How well can we recover the angular power spectrum of the cosmic 
microwave background (CMB) and the cosmological parameters from maps of a given 
experiment, which can have holes and regions with no n- uniform sampling? How well 
can we subtract foregrounds from a multifrequency map? An important issue for 
small angular scale simulations is their great computational cost. It is important to 
have fast and accurate methods to do simulations. This paper is concerned with the 
comparison of two of these methods. 

Large-scale simulations of the microwave sky are based on the expansion of the 
temperature contrast in spherical harmonics, 

Xrp lmax m=+l 

y-(M)= E E a lm Y lm (0,<P)- (1) 

1=1 m=—l 

Typically l max < 40 in COBE-like maps, since the number of coefficients, a/ m , to 
be calculated is (l ma x + I) 2 ; an d there are questions of numerical accuracy in high 
order calculations. In order to make high resolution maps — resolutions ~ 10'— 
l max must be of order 10 3 so that the number of coefficients to be calculated is of 
order 10 6 ; moreover, the 10 6 spherical harmonics must be evaluated accurately in 
roughly 10 6 locations in a full-sky map (10 4 locations in a 20° x 20° field). Thus, 
for high angular resolution, the spherical harmonic expansion becomes difficult even 
with modern computing power. Although simulations based on spherical harmonics 
are feasible (Hinshaw, Bennett & Kogut 1995; Kogut, Hinshaw & Bennett 1995; 
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Jungman et al. 1995), the Fast Fourier Transform (FFT) is shown to be useful in 
order to make simulations efficiently on small angular scales. Hereafter, simulations 
using spherical harmonics are called "SH simulations", while those based on plane 
waves are called "FFT simulations". The main goal of this paper is the comparison 
of the accuracies and the CPU costs corresponding to SH and FFT simulations. Both 
approaches use the C e = J2m=-e \ a em\ 2 /(^+ 1) quantities corresponding to the model 
described in Sec. 2. The method used to extract the angular power spectra from 
the simulated maps is described in Sec. 3. FFT simulations are described in Sec. 4, 
while SH simulations are presented in Sec. 5. Comparisons between FFT and SH 
simulations are given in Sec. 6 and, finally, Sec. 7 contains the main results and a 
general discussion. 

2 The model 

All the simulations presented in this paper are based on a single cold dark matter 
(CDM) model for large scale structure formation. In this simple CDM model, no 
reionization modifies the CMB anisotropies, the background is flat (Qo = 1), the 
cosmological constant vanishes, scalar fluctuations are Gaussian, their spectrum is 
scale- invariant, and tensor fluctuations are absent. This model is hereafter called the 
minimal CDM model. The free parameters of this model are: the amplitude of the 
angular power spectrum of the CMB, the density parameter of the baryonic matter 
Q B and the reduced Hubble constant h. The parameters Q B and h are involved in 
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the transfer function; their values are assumed to be 0.03 and 1/2, respectively. 

In the above model, the normalization of the angular power spectrum is fixed by 
to rms quadrupole of the CMB anisotropy generated by scalar modes by using the 
estimator Q rm s-ps found by fitting the observed temperature fluctuations assuming 
a scale-invariant primordial density power spectrum. In the absence of tensor modes, 
experiments measuring at large angular scales, such as COBE (Smoot et al. 1992, 
Bennett et al. 1992, Wright et al. 1992) and TENERIFE (Hancock et al. 1994), 
provide estimations of Q rm s-ps- The Ce coefficients have been taken from Sugiyama 
(1995) and renormalized according to the four- year COBE data (Q rm s~ps — 18 [iK, 
Gorski et al. 1996). All the simulations performed in this paper are based on these 
Cn coefficients. 



3 Obtaining the angular power spectrum from sim- 
ulated maps 

In this section, we describe the method used in order to extract the angular power 
spectrum from a given map. Let us give some useful definitions. 

The autocorrelation function can be defined as C(9) = C a=0 (9), where 

The angle between the unit vectors Hi and n 2 is 9. The angular brackets stand 
for the mean over many full realizations of the CMB sky. The quantity (^) (7 (n) is 

3 



the temperature contrast in the direction n after smoothing with a Gaussian beam 
described by a = 0A259 FWHM . 

Function C a (9) can be expanded in the following form: 

1 oo 

C M = ^E( 2 ^+ l)C e Pe(cos6)e^+ ^ 2 . (3) 

47r 1=2 

where the Cg quantities define the power spectrum of the CMB. These quantities have 
been calculated in many theoretical models of structure formation. 
From Eq. (3) one easily obtains the relation: 

C e (a) = e^ +0 - 5 ^ 2 C e = ^ j o C a (9)P e (cos9) mnOM . (4) 

Given a simulated map, the function C a {9) is obtained by using Eq. (2). For each 
value of 9, many pairs of direction vectors (ni, n 2 ) — forming a fixed angle 6 — are 
randomly chosen. To be more precise, the direction vector n\ is chosen to point at a 
random node of the map — where the temperature is known — while the direction n 2 is 
randomly placed on a cone centered on the node corresponding to n\. The half-width 
of this cone is 9. The temperature corresponding to the direction n 2 is obtained by 
interpolating between neighboring nodes. As discussed below, mathematical interpo- 
lation introduces only a small error in the resulting spectrum for values of 9 lying in 
the interval (9 m i n , 9 max ), where 9 m i n is the angle separating two neighboring nodes 
and 9 max is an angle to be experimentally obtained. This angle must be: (a) smaller 
than the size of the simulated region in order to allow us to place a great number of 
pairs (fti, n 2 ) inside it, and (b) large enough to include as many scales as possible. 
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Since we know the values of the function C a {9) only interval (9 m i n , 9 max ), we 
cannot compute C^(cr). Instead, we have calculated the quantities 

D e {a) = J0^ 2 C a {9)P,{cos9) sm9d9 . (5) 

These quantities are to be compared with the corresponding theoretical ones, which 
can be obtained as follows: First, from Sugiyama's Cg coefficients of the minimal 
CDM model and from Eq. (3) (with the sum extended from i = 2 to £ = 1100), the 
autocorrelation function is estimated in the interval (9 min , 9 max ). This theoretical 
autocorrelation function plus Eq. (5) then allows us to compute the theoretical values 
of D e (a). 

It is hereafter said that the quantities D e (a) define the modified power spectrum, 
which is different from the true angular power spectrum given by the coefficients 
Ci{<t). 

In all the simulations we are considering in this paper, the following values have 
been selected: 9 min = 5', 9 max = 4.5° and 9 FWHM = 10'. These values are similar to 
those of modern projects for the detection of small angular scale anisotropics. 

A very accurate estimation of the mean involved in Eq. (2) requires multiple 
realizations of the microwave sky. In the case of a single full realization, a certain 
error in this mean — producing an error in C a {6) — is unavoidable. This error is usually 
measured by the cosmic sample variance. The smaller the sky coverage, the greater 
the errors in the resulting autocorrelation functions and spectra (Scott, Srednicki 
& White, 1994). As shown below, the spectra resulting from a small 20° x 20° 
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coverage are rather different than the theoretical ones. Although coverages greater 
than 20° x 20° cannot be achieved by using FFT simulations, we can use multiple 
disconnected 20° x 20° regions in order to improve on the estimate of the mean involved 
in Eq. (2), for values of 6 smaller than 9 max . In other words, FFT simulations become 
useful as a result of the fact that various independent small coverages of the sky can 
be used in order to obtain a good estimate of the autocorrelation function for small 
angular scales; however, neither a 20° x 20° map nor an ensemble of them has good 
information about angular correlations larger than 9 max . 



4 Simulations based on the FFT 

Our FFT simulations are based on the use of the standard Ce coefficients, which have 
already been calculated for many cases. As pointed out by Bond and Efstathiou 
(1987), a FFT simulation is only possible in the case of anisotropies corresponding 
to large ^-values and small regions of the sky. A quantitative study of the validity 
and usefulness of this type of simulation has not previously been presented. In order 
to do this study, we compare results from FFT simulations with results from SH 
simulations. 

The temperature contrast in the simulated region can be obtained from the fol- 
lowing formula: 

ST N 

^-(M)= E D(e 1 ,e 2 )e- i ^ + ^ (6) 

S1,S2 = -N 

where t\ = 2-ksi/ A and i 2 = 2irs 2 /A, A being the angular size of the elemental square. 
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The distribution of the quantities D(£\, £ 2 ) is assumed to be Gaussian, with zero mean 
and with variance C^/A 2 , where £ = (£\ + £\) x ^ 2 ■ The quantities D(£ 1 ,£ 2 ) satisfy the 
relation D(-£ u -l 2 ) = D*(£ 1 , £ 2 ). 

The use of plane waves requires negligible spatial curvature in the simulated re- 
gions. This restricts our simulations to small parts of the microwave sky. In this 
paper, it is verified that a region of 20° x 20° (~ 1% of the total sky) can be success- 
fully mapped by using the FFT. The value 9 max = 4.5° is marginally compatible with 
the size of this region (many pairs (Hi, H 2 ) forming an angle of 4.5° can be located 
inside it). An angular scale of 4.5° corresponds to £ = 40 and, consequently, the 
spectra obtained from 20° x 20° FFT simulations is only estimated for £ > 40. In 
practice, the accuracy in the estimation of the power spectrum from these simulations 
seems to be good for £ > 100 (see Sec. 6). 

In a 20° x 20° region, the assumption of flatness produces a deformation of the true 
maps, but this deformation is not expected to hide their main features. We remind 
the reader that only correlations at scales smaller than 4.5° are being estimated and, 
on these scales, the curvature is clearly negligible. 

Since the FT involves periodic boundary conditions, our simulations are not phys- 
ically significant at all points of the elemental square. Only the simulation of the 
central region is physically admissible. On account of this fact, Fourier transforms 
are performed in a large ~ 40° x 40° square (A = 40°) and the resulting simulation 
is truncated to the central ~ 20° x 20° region. The curvature could be important 
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if we wanted to get a physically significant mapping of the large square; however, 
the large square is only an auxiliary element without any physical significance. The 
number of points on each edge of the large square is taken to be 2N=512. This choice 
corresponds to 9 min ~ 5'. Simulations using an ~ 80° x 80° square with 1024 points 
per edge yield similar results; hence, the use of more than 512 points per edge is not 
necessary. 

On an IBM 30-9021 VF, the CPU cost of our FORTRAN code is ~ 5.5 minutes 
per simulation. 

Figure (1) shows a FFT simulation of a ~ 20° x 20° region. No unusual features 
are apparent near the edges. The numbers appearing in this Figure are the minimum 
and maximum values of the temperature contrast in our simulation. Only a few points 
of the simulation approach these values. 

5 Simulations based on spherical harmonics 

Our SH simulations are based on Eq. (1) with £ max = 1100 and £ min = 40. The ae m 
coefficients have been generated as statistically independent random numbers with 
variance (|c^ m | 2 ) = Cee^^ - 5 ^" 2 and zero mean. The spherical harmonics have been 
carefully calculated. These simulations include scales smaller than those considered 
in previous simulations (see Hinshaw, Bennett & Kogut 1995; Kogut, Hinshaw & 
Bennett 1995; Jungman et al. 1995). The small values of 9 min and a considered in 
our simulations require the use of large values of £ giving information about small 
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angular scales. 

In this paper, SH simulations are performed only for comparison with FFT simu- 
lations; for this reason, some features of our numerical SH simulations are identical to 
those of the FFT simulations of Sec. 4. In particular, our SH simulations are extended 
to 20° x 20° regions of the sky. These regions are assumed to be uniformly covered 
and the number of points per edge is 256 (9 min ~ 5'). A smoothing corresponding to 
FWHM = 10' has been performed, and the angle 9 max = 4.5° has been assumed. 

Simulations of the two kinds defined up to now have been obtained for the minimal 
CDM model. The Ce coefficients have been obtained and normalized as detailed in 
Sec. 2. 

On an IBM 30-9021 VF, the CPU cost of our FORTRAN code is ~ 53 minutes 
per simulation. This cost is about ten times greater than that of the FFT simulations 
described in Sec. 4. 

SH simulations of ~ 20° x 20° regions look like the FFT simulation of Fig. 1. 
Differences between SH and FFT simulations can be noticed only by comparing the 
spectra obtained from them. These comparisons are discussed in the next section. 

6 Results 

Let us now compare the theoretical Dg(a) quantities with those obtained from FFT 
and SH simulations. 

In all the panels of Fig. (2), the solid lines correspond to the modified theoretical 
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spectrum and the dashed lines show the modified spectrum obtained from simulations. 
The left (right) panels of this figure show the modified spectra extracted from FFT 
(SH) simulations. The top, intermediate, and bottom panels exhibit the modified 
spectra obtained from 1, 13, and 52 small simulations of 20° x 20° regions, respectively. 
The total area covered by the 1, 13, and 52 small simulations is 400, 5400, and 21600 
( ~ \ of the full sky area) squared degrees, respectively. Hereafter, N s stands for the 
number of 20° x 20° simulations. 

In order to measure the deviations between the theoretical and simulated spectra 
in each case, the following quantities have been calculated and presented in Table 1 : 
The mean, Ml, of the quantities 0.69£(£ + l)D e (a) x 10 10 (column 3); the mean, M2, 
of the differences between the theoretical and simulated values of these quantities 
(column 4); the typical deviation, S, of the same (column 5); and the mean, MA, of 
the absolute value of these differences (column 6). These quantities are estimated in 
appropriate £ intervals (column 7). 

As can be seen from the top panel of Fig. 2 and from Entries 1 and 2 of Table 1 
(see MA and E), a single 20° x 20° map leads to a poor estimation of the angular 
spectrum in both FFT and SH simulations. A greater sky coverage is necessary. The 
intermediate panels of Fig. 2 and Entries 3 and 4 of Table 1 prove that the use of 
thirteen 20° x 20° maps strongly improves on the resulting spectra (see MA and S). 
The quantities MA and £ decrease as N s increases. The total decrease undergone 
by MA and E as N s goes from 1 to 13 (see the top and intermediate panels of Fig. 
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2 and Entries 1-4 of Table 1) is larger than that appearing as N s increases from 13 
to 52 (see the intermediate and bottom panels of Fig. 2 and Entries 3-6 of Table 1). 
The quantity M2 is of order 10~ 2 whatever the value of N s may be. The absolute 
value of M2 is always much smaller than that of MA; this means that the deviations 
between the theoretical and the simulated spectra are essentially oscillatory. This 
fact strongly supports the idea that moderate coverages could give good spectra after 
removing oscillations. The above results can be independently obtained from both 
SH and FFT simulations. 

For N s = 52 and the ^-interval (40,100), the quantities \M2\, MA, and £ cor- 
responding to FFT simulations are greater than those of SH simulations (compare 
Entries 7 and 8 of Table 1). This means that, for £ < 100, SH simulations are 
preferable (in spite of their large CPU cost). 

The above results have been verified in other realizations of FFT and SH simula- 
tions corresponding to the N s values 1, 13 and 52; however, only one set of realizations 
is presented (Fig. 2) for the sake of brevity. 

Although this paper is devoted primarily to the comparison of FFT and SH simu- 
lations in the absence of noise, the effect of uncorrelated noise at the level of 27.3 fiK 
in a 20° x 20° region has been estimated. In order to do this, a 20° x 20° sim- 
ulation involving only uncorrelated noise has been made. The spectrum has been 
obtained from the resulting map following the prescriptions of Sec. 3. The quanti- 
ties 0.69£(£+ i)De(a) x 10 10 appear to be smaller than 10~ 2 for all I. These values 
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are negligible compared to those corresponding to the fluctuations of the microwave 
background (see Fig. (2); 0.69£{£+ 1) D e (a) x 10 10 ~ 1). This proves that the effect of 
uncorrelated noise in our 20° x 20° maps can be neglected. This fact has been verified 
by superimposing the above white noise map upon a SH map. The spectra obtained 
from the resulting map is practically the same as in the absence of noise. 

7 Conclusions and discussion 

Small scale simulations of the microwave background are expected to be very useful 
in order to analyze observational data from both current and future experiments. For 
scales smaller (greater) than ~ 1.8° (£ ~ 100), the accuracy of the FFT simulations 
appears to be comparable to (worse than) that of the SH simulations. The CPU cost 
of the SH simulations is about ten times greater than that of the FFT simulations. 
These conclusions strongly support the use of FFT simulations as a fast and accurate 
tool leading to very good spectra in the case of small angular scales (£ > 100). 

There are various sources of the deviations between the theoretical spectrum and 
the spectrum extracted from FFT simulations. Significant contributions to these 
deviations are expected from partial coverage, boundary conditions in FFT maps, and 
the use of mathematical interpolation in order to assign a temperature contrast to the 
direction n 2 (see Sec. 3). Such interpolation introduces non-physical information even 
for the most sophisticated interpolation methods. An error in the resulting spectrum 
associated with this fact seems to be unavoidable. Fortunately, this error appears 
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to be small from 9 m i n to 6 max . Boundary effects associated with the FFT have been 
minimized by simulating the core of an auxiliary 40° x 40° region. 

Uncorrelated noise does not appear to be relevant for the coverages considered in 
this paper. 

Since we are only considering angular scales smaller that 4.5°, curvature effects 
are expected to be negligible. 

The main conclusions of this paper are independent of admissible renormalizations 
of the Ci coefficients. 
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Figure Captions 

Fig. 1. FFT simulation of a 20° x 20° region of the microwave sky. Numbers 
correspond to the minimum and maximum values of the temperature contrast. 

Fig. 2. Each panel shows the quantity O.Q9£(£ + l)De(a) x 10 10 as a function of 
log(£) in various cases. In all the panels, the solid line shows the theoretical modified 
power spectrum and the dashed line shows the modified spectrum extracted from 
simulations. All the modified spectra have been obtained for a = 10', 6 min = 5', and 
Qmax = 4.5°. Top, intermediate, and bottom panels correspond to the N s values 1, 13, 
and 52, respectively. Left (right) panels show modified spectra obtained from FFT 
(SH) simulations. 
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TABLE 1 



COMPARING THEORETICAL AND SIMULATED SPECTRA 



Type of 
simulations 


N s 


Ml 


M2 




E 




MA 




^-interval 


FFT 


1 


1.07 


-1.01 x 10~ 


-2 


7.13 x 10" 


-3 


1.58 x 10" 


-i 


40- 


1000 


SH 


1 


1.07 


-1.33 x 10~ 


-2 


6.96 x 10" 


-3 


1.76 x 10- 


-i 


40- 


1000 


FFT 


13 


1.07 


-1.13 x 10" 


-2 


2.58 x 10" 


-3 


5.90 x 10" 


-2 


40- 


1000 


SH 


13 


1.07 


-1.85 x 10- 


-2 


3.08 x 10- 


-3 


8.00 x 10- 


-2 


40- 


1000 


FFT 


52 


1.07 


-2.04 x 10- 


-2 


2.25 x 10- 


-3 


5.21 x 10- 


-2 


40- 


1000 


SH 


52 


1.07 


-2.56 x 10- 


-2 


2.30 x 10- 


-3 


6.06 x 10- 


-2 


40- 


1000 


FFT 


52 


1.24 


-5.29 x 10- 


-2 


2.46 x 10- 


-2 


1.61 x 10- 


-1 


40 - 


- 100 


SH 


52 


1.24 


-1.54 x 10- 


-2 


1.05 x 10" 


-2 


6.54 x 10" 


-2 


40 - 


- 100 
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